function [x, acc] = gaussseidel(A,b,n,x0)

% [x, acc] = gaussseidel(A,x0,b,n)
%
%   Solves Ax = b using the Gauss-Seidel method.
%
%   x0 := initial x
%   n  := # of iterations

L = tril(A);
Linv = inv(L);
U = triu(A, 1);

G = Linv * U;
c = Linv * b;
x = x0;

acc = zeros(3, n);

for k = 1:n
    
    x = c - G*x;
    acc(1, k) = k;
    acc(2, k) = x(1);
    acc(3, k) = x(2);
    
end